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. Aims. The long term magnetohydrodynamic stability of magnetized transonic/supersonic jets is numerically investigated using a spatial 
^ approach. We focus on two-dimensional linearly-unstable slab configurations where the jet is embedded in a flow-aligned uniform magnetic 
1 1 1 field of weak amplitude. We compare our results with previous studies using a temporal approach where longitudinally periodic domains were 
adopted. 

Methods. The finite-volume based versatile advection code is used to solve the full set of ideal compressible MHD equations. We follow the 
development of Kelvin-Helmholtz modes that are driven by a white noise perturbation continuously introduced at the jet inlet. 
Results. No noticeable difference is observed in spatial simulations versus analogous temporal ones during the linear and early non-linear 
evolution of the configuration. However, in the case of transonic flows, a different long-term scenario occurs in our spatial runs. Indeed, after 
the large-scale disruption of the flow, a sheath region of enhanced magnetic field encompassing the jet core forms along the whole flow. This 
provides a partial stabilization mechanism leading to enhanced stability for later times, which is almost independent of the initial magnitude of 
the magnetic field. The implication of this mechanism for the stability of astrophysical jets is discussed. 
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;_i 1 . Motivation systems is now well established from observations (Ray et al. 

1997; Pushkarev et al. 2005; Gabuzda et al. 2004) and jet- 

Astrophysical jets are observed to propagate over very long dis- launching models (Casse & Keppens 2Q02j 2004 and refer . 

tances with respect to their radial extent, typically up to one ences therein) The magnetic field is probably too weak to pre- 

thousand times the jet radius. This is, for example, the case of yent Qr strongly weaken ±e Hnear development of the KH in _ 

many jets emanating from young stellar objects (YSO) and ac- stabmty (Baty 2005; Appl & Camenzind m2) . However, the 

live galactic nuclei (AGN). These observations seem to be in evolution of the MHD kh instability allows a much richer 

disagreement with the magnetohydrodynamic (MHD) stability complexity in the non _ii near phase, compared to its purely hy- 

theory and numerical simulations, which predict the strong de- drodynamic counte rpart. For example, a weakly magnetized 

velopment of instabilities that consequently threaten the jet's shear flow configuration has three dynarmca n y diffe rent sub- 

colhmation. Particularly important is the Kelvin-Helmholtz regimes depending on the value of the Alfven Mach number 

(KH) instability, typical of shear flows, which is shown to dis- (Raty et al 2003; Jonefj et al ml) A recent numerical in . 

rupt the jet ona few sound transit timescales (see reviews by vestigation showed that the overall scen ario, even if it is en- 

Birkinshaw[199ll and Ferrari 1998). This is mostly evident in riched by magnetic reC onnection events and large-scale coales- 

purely hydrodynamic simulations, where supersonic flows un- cence eflfectSi leads tQ quick dismption in MHD as well (Baty 

dergo a strong turbulent disruption driven by the KH modes & Keppens 2006i referred to as Paper l in this study) The latter 

(Bodo et al. 1995, 1998). result is valid for the non-linear evolution of uniformly magne- 

Among the possible mechanisms that are likely to have tized j ets having wide i y different plasma-beta, sonic/Alfvenic 

a stabilizing effect is the role of magnetic fields. Indeed, the Mach number s, and shape of the flow profile. However, a tem- 

presence of a large-scale magnetic field in accretion-disk-jet poral approac h with periodic boundary conditions imposed at 

the inflow/outflow boundaries is used in Paper I. Indeed, this 
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temporal procedure is preferable when one wants to focus on 
long-term temporal evolution, as it allows one to compute the 
time evolution of a jet section at relatively low numerical cost. 
On the other hand, the temporal approach has several limita- 
tions (Baty et al.|2003). The most important one is that the in- 
teraction between different parts of the jet is not well taken into 
account, as the convective effect of the mean flow is ignored in 
periodic configurations (even in cases of extended domains). 

Our aim is to investigate the development of the KH insta- 
bility using a spatial approach, where we follow the evolution 
of a perturbation introduced at the jet inlet of a non-periodic 
configuration. This implies the use of an extended numerical 
domain, as the study is limited by the transit time through the 
jet configuration. We adopt similar physical parameters for the 
jet as those assumed in Paper I, in order to make a close com- 
parison with results obtained using a temporal approach. We 
consider only two dimensional (2D) configurations, mainly be- 
cause the numerical cost of a similar study in three dimensions 
(3D) is beyond our current computational capabilities. 

We mainly focus here on weakly magnetized transonic 
flows. This is a relevant regime for studying the early phase of 
jet propagation. Indeed, before reaching the asymptotic region 
(where the flow is believed to be supersonic), the jet is acceler- 
ated in an intermediate region where it must survive instabili- 
ties. And the latter intermediate regime is of prime importance 
for KH modes, as the spatial growth rates are larger for tran- 
sonic flows vs supersonic flows. A selected supersonic jet case 
is nevertheless investigated to address the generality /limitation 
of our transonic results. Previous numerical studies have al- 
ready adopted a spatial approach (e.g. Hardee et al. 1992; 
Hardee et al. 1995; Micono et al. 119981 1. However, the latter 
studies focus mainly on relatively strongly magnetized jets and 
describe the whole jet in order to make comparisons with ob- 
servational characteristics; and for this purpose they necessar- 
ily sacrifice resolution in the description of modes (especially 
in 3D). Finally, note that we do not deal here with magnetic 
current-driven instabilities, which are driven by the presence 
on an electrical current density parallel to the magnetic field in 
3D helical field geometry. 

This paper is organized as follows. We present the jet model 
and numerical setup used in this study in Sect. 2. In Sect. 3, we 
show the numerical results, which are split into an overview of 
a typical transonic flow run, a comparison with temporal sim- 
ulations, an in-depth study of a partial stabilization mechanism 
observed in our simulations, and an overview of a supersonic 
run. Finally, conclusions are drawn in Sect. 4. 

2. Jet model and numerical setup 

In the present study, we use the full set of ideal compress- 
ible MHD equations, as we consider perfectly conducting plas- 
mas having very large kinetic and magnetic Reynolds numbers. 
This description is a suitable proxy for astrophysical jets envi- 
ronments. Thus, we rely on the inherent numerical resistivity 
and viscosity (due to scheme discretizations) to mimic the dis- 
sipative processes. It is generally admitted that this approach 
provides an adequate model for subgrid-scale dissipation (see 
Paper I, and the discussion in Jones et al. 1997). However, a 



Table 1. Numerical and physical parameters of our simulations. 
Run 2b corresponds to the numerical linear convergence study 
and is done on a smaller domain. 



Run 


M s 


M A 


M f 


Domain 


Resolution 


1 


1 


5 


0.98 


[0,40]x[-8,8] 


1200x480 


2 


1 


7 


0.99 


[0,40]x[-8,8] 


1200x480 


2b 


1 


7 


0.99 


[0,20]x[-2,2] 


200x40, 300x60, 












600x120, 800x160 












1200x240,1600x320 


3 


1 


10 


0.99 


[0,40]x[-8,8] 


1200x480 


4 


1 


14 


0.997 


[0,40]x[-8,8] 


1200x480 


5 


1 


5 


0.98 


[0,80]x[-16,16] 


1200x480 


6 


1 


7 


0.99 


[0,80]x[-16,16] 


1000x600 


7 


1 


14 


0.997 


[0,80]x[-20,20] 


1000x600 


8 


3 


7 


2.75 


[0,80]x[-16,16] 


1000x600 



The numerical domain is rectangular, of size [0, L x ] x [-L Y /2, L y /2]. 
In transonic runs, the amplitude of the jet velocity is V = 1-29 and the 
initial magnetic field intensity is varied. In the supersonic run, Vq = 
3.87. 



convergence study by repeating the simulations using different 
grid resolutions is recommended. 

The set of MHD equations is solved in 2D Cartesian geom- 
etry, with the direction of flow propagation being along the x 
axis. The transverse coordinate is y. We consider a jet having 
an initial background flow profile: 
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where Vo is the velocity amplitude of the jet. The parameters 
Rj and F control the shape of the velocity profile; Rj is the jet 
radius, and F determines the thickness e of the shear flow (e = 
Rj/F, see Paper I). For convenience, we also introduce the half- 
thickness of the vorticity layer a = e/2. In this study, we use 
a flow profile similar to the Pr2 case studied in Paper I ; more 
explicitly we take Rj = 0.5, and F = 2 (thus a = 0.125 in our 
case, which differs from Paper I where a = 0.05 for all profiles). 
We consider a medium with initially uniform temperature To = 
1 and density po = 1, thus defining our normalization; the sonic 
speed c s = (yPo/po) 1 ^ 2 is then equal to 1.29 in our units (with 
the adiabatic index y = 5/3). Consequently, the time needed for 
sound waves to travel through the jet radius is Rj/c s = 0.3876 
in our units. 

The flow is embedded in an initially uniform longitudinal 
magnetic field of amplitude Bq. The parameters used in the dif- 
ferent runs are listed in Table [T] where we take the following 
definitions for the sonic and Alfvenic Mach numbers, M s = 
Vq/c s and M A = Vq/V a (V a = Bq/ -y^oo is the Alfven velocity). 

The fast Mach number is defined as Mf = Vq/ Jcj + V^.As 
already emphasized, we focus on transonic cases (M s = 1). 
However, a supersonic case M s — 3 is also considered for com- 
parison and for the sake of completeness in Sect. |3.4| Relatively 
weak magnetic field amplitudes are taken with M A in the range 
[5-14] which correspond to values of the so-called "weak field 
regime" or "disruptive regime" (e.g. Jones et al. 1 19971 Baty et 
al. 2003). Indeed, higher M A values seem to be unrealistic for 
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Fig. 1. Density distribution in the x-y plane obtained at time t — 16 for run 2 (see Table [TJ. Dark regions correspond to low 
density values. A linear scale is used with density values ranging from 0.6 to 1.2. The unit of length is the jet diameter. 
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Fig. 2. Same as Fig.Q] at t = 22. 
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Fig. 3. Same as Fig.Q] at t = 34. 



astrophysical jets and lower values lead to stable configurations 
(also probably unrealistic). 

We use the general finite-volume based versatile advection 
code VAC (T6th [T996i l for our simulations. The MHD equa- 
tions are solved by selecting the explicit one-step total varia- 
tion diminishing (TVD) scheme with minmod limiting (Colella 
& Woodward 1984; Harten 1983). This is a second-order ac- 



curate shock-capturing method using a Roe-type approximate 
Riemann solver. Our VAC simulations also apply a projection 
scheme at every time-step to remove any numerically gener- 
ated divergence of the magnetic field up to a predefined ac- 
curacy. Spatial simulations require an extended numerical do- 
main; Table Q] summarizes the different domain sizes used. In 
short, we employed two different length values for the numer- 
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ical domain: a "short" one with L x = 40 having a relatively 
high spatial resolution, and a "large" one with L x = 80 having 
a lower resolution to study the large-scale stability of the jet. 
As an initial condition, the jet profile (1) was set up across the 
whole computational domain. The jet continuously enters the 
domain by its left boundary and free outflow conditions were 
used on all other boundaries (by imposing a zero gradient for 
every variable). Lateral boundaries were taken far enough away 
from the jet flow to avoid spurious effects. Thus, at the jet in- 
let (left boundary at x — 0), the background velocity profile 
(Q]i is maintained in time. Moreover, to excite the unstable KH 
modes, we continuously perturbed the jet inlet (at x = 0) by 
adding a small amplitude transverse velocity, which is taken 
to be one percent of the magnitude of the background flow. 
Specifically, as Paper I, we imposed a white noise perturbation 
by using a random number generator. This is a valid proce- 
dure for selecting the linearly fastest growing mode as shown 
by Zhao et al. (1992). Concerning the other physical variables 
at this boundary, we imposed fixed conditions. Indeed, the real 
physical conditions are not known and this latter, somewhat ar- 
bitrary, choice is convenient for our numerical procedure. 

3. Results 

Before examining in detail the simulation results, it is informa- 
tive to have an overview of the typical evolution of a transonic 
case (run 2 in Table 1). 

3.1. Overview of a transonic run 

Clearly, two simulations will always differ in their details, in 
particular due to the random character introduced by the pertur- 
bation. However, the following scenario can be inferred from 
our simulations. Figures [HE] and [3] show the first stages of the 
transonic run 2 (see Table 1). First, one can clearly see in the 
density snapshots of Fig. 1 the initial development of a few un- 
stable KH wavelengths, which take the form of two rows of 
vortices located on the two shear layers of the mean flow. We 
term this set of growing instabilities that are convected along 
the flow direction, the instability train. Later, these vortical 
structures undergo pairing-merging events (due to an inverse 
cascade effect towards large scales, as described in Paper I) ac- 
companied by magnetic reconnection, and finally a strong dis- 
ruption (Figs. ISO. This scenario is characteristic of the disrup- 
tive regime, and is expected from temporal studies. In the next 
section, a detailed comparison is made with results obtained in 
a similar temporal simulation, to quantify the convective effect 
in our spatial run. 

The disruption occurring at t = 34 in Fig. [3] completely de- 
stroys the initial laminar structure of the jet's flow. This is il- 
lustrated in Fig. [4] where the normalized x-momentum profile 
at the location of the large-scale disruption in Fig.[3](->c ~ 20 at 
t = 34) and at the inlet is plotted. As a consequence, we can es- 
timate a maximum propagation distance over which the jet flow 
is able to maintain its coherence, that is roughly equal to 40 jet 
radii in our simulation. This result is consistent with a value 
deduced from previous temporal studies, where a similar value 
for a disruption length was obtained by translating the typical 




j j 



y 

Fig. 4. One dimensional cut of the normalized x-momentum 
(momentum pV divided by the initial momentum po Vq) profile, 
obtained at x = 20 (dashed line) and x = (plain line) in run 2 
at t = 34 (Fig.D. 



disruption time (Paper I). This is also in close agreement with 
the value deduced from the additional temporal run performed 
in the present study (see Sect. |3.2| >. 

In Figs. [2] and [3] one can see two additional striking fea- 
tures. First, at the head of the instability train, one can see the 
formation of a propagation front (located at x » 24 in Fig. [2] 
and at x ~ 37 in Fig. [3j, which is found to travel in the down- 
stream direction at approximately the jet's speed. Behind this 
front, a sinusoidal density pattern is also formed (Fig. [3] be- 
tween i * 23 and x * 37). The latter structure arises from 
secondary instabilities which are driven at the head of the in- 
stability train. This "convective generation" of instabilities is 
discussed in Yamamoto et al. 1988. Indeed, the primary modes 
of the instability train can non-linearly trigger secondary KH 
modes when convected by the mean flow, since the downstream 
medium is taken to be linearly unstable. With our jet-like ve- 
locity profile, this leads to the formation of a single vortex-like 
street (sinusoidal pattern) located in the jet core, as one can see 
in the density snapshot of Fig 3. We verified that the latter phe- 
nomenon has no effect on the large-scale disruption of the jet 
(occurring far behind the front), and thus no consequence on 
the flow survival. The second striking feature is the formation 
on the two edge layers of a trail structure, which is found to 
propagate in the upstream direction toward the jet's inlet (Fig. 
[3]between x ~ 7 and i s 11). This trail is also visible at early 
times, between i«8 and x ~ 10 in Fig.|2](but less clearly). We 
showed that this particular phenomenon is associated with a lo- 
cal magnetic amplification process (during the formation of the 
KH vortices of the instability train) which is able to reach the 
jet inlet at later times. A partial stabilization mechanism thus 
ensues as demonstrated in detail in Sect. 13.31 
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3.2. Comparison with the temporal approach 

We compare here the results obtained in our typical transonic 
spatial case (run 2) with results obtained in a similar temporal 
simulation. Note that the physical parameters of the jet con- 
figuration are the same, and also that we use the same code, 
VAC. The only difference comes from the numerical procedure, 
where periodic boundary conditions are assumed in the longi- 
tudinal direction for the temporal run (see Paper I for more de- 
tails). The periodic length of the domain is taken to be L x = 8, 
to allow a few wavelengths of the linearly fastest growing mode 
to fit the extent of the spatial domain (see Fig. [6}. 

Before examining the results in details, one needs to have 
in mind the main mathematical differences between the two ap- 
proaches (spatial and temporal) for the description of linearly 
growing instabilities in such flowing plasmas. The standard lin- 
ear analysis assumes an infinitely long jet in equilibrium. The 
linear perturbations are thus conveniently expanded using the 
classical normal mode form: 



f(x,y, t) - A(y) exp U(kx - wf)], 



(2) 



for any perturbed physical quantity /. Here, u> and k are respec- 
tively the frequency and wavenumber of a given normal mode. 
The linearized MHD equations then yield a dispersion relation: 



D(k, to) = 0, 



(3) 



which is specific to the jet configuration. The stability of nor- 
mal modes is obtained by solving this dispersion equation. This 
can be done in two ways. First, in a temporal approach, a com- 
plex frequency u - o> r +iy and real wavenumber k are assumed. 
Consequently, y is the temporal linear growth rate of the corre- 
sponding wavenumber k. On the other hand, a spatial approach 
considers that the frequency to is real and that the wavenum- 
ber k = k r + zT is now complex. As a consequence, the linear 
spatial growth rate is directly obtained by F, and r gives the 
corresponding characteristic growth length. This latter formal- 
ism better accounts for the convective nature of the KH insta- 
bilities, which are not growing at a fixed point in space, but 
instead are convected along the fluid forming a spatially grow- 
ing pattern. Both approaches can be connected using the group 
velocity v g of the modes, with T = y/v g (Drazin & Reid 1981; 
Appl & Camenzind 1992). Moreover, as the modes under con- 
sideration in this study are weakly dispersive, one can also use 
the relation T y/v p where v p is the phase velocity (Payne & 
Cohn [T985l Appl 1996). 

Figure|5]shows a zoom of the instability train (plasma den- 
sity) obtained in the typical transonic spatial simulation dis- 
played in Fig. Q] One can see that the latter snapshot is very 
similar to the density map obtained from a temporal simulation 
(Fig. [6]). In both runs, a double layer of three vortices localized 
on the shear flow is clearly visible. The longitudinal extent of 
each vortex is A « 2 in our units, which corresponds to a nor- 
malized wavenumber ka w 0.4. This is in excellent agreement 
with results obtained from a linear stability analysis of the slab 
jet where the wavenumber of the linearly fastest growing mode 
is ka * 0.4, and does not depend significantly on the form fac- 
tor F of the jet (e.g. the growth rates curve in Fig. 1 of Paper I, 




Fig. 5. Zoom on the instabilities train, taken from the spatial run 
2 at t — 16 (Fig. |TJ. A linear scale is used with density values 
ranging from 0.6 to 1.2. The unit of length is the jet diameter. 




Fig. 6. Density distribution for a temporal simulation obtained 
at a similar time as Fig.|5](see text). The unit of length is the jet 
diameter. 



where the maximum rate occurs for k « 8 with a profile having 
a = 0.05). The wavelength of the fastest growing mode of the 
slab jet is also very close to the wavenumber obtained for a sin- 
gle shear layer embedded in a uniform magnetic field, where 
ka * 0.41 (Keppens et al. 1999; Miura & Pritchett 1982). This 
is not surprising as the aspect ratio of the jet flow profile con- 
sidered in this study is quite large (Rj/a = 4 » 1), leading thus 
to a rather weak interaction between both layers during the lin- 
ear regime. A similar conclusion was previously reported by 
Min(1997). 

In order to have a quantitative comparison of both ap- 
proaches, we measured the growth rates by following the 
growth of the mean transverse kinetic energy: 



Eky(x, t) 



— I < 

Ly J-Ly/2 



1 2 

dy-p(x,y,t)v y (x,y,t), 



(4) 
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Fig. 7. Isocontours of the mean transverse kinetic energy (see 
Eq. (01) in the (t,x) plane. Ten linearly spaced isocontour val- 
ues (between a minimum and a maximum) are considered. The 
dashed line is the characteristic curve corresponding to an in- 
stability starting to grow at t — and x — 0. 
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Fig. 8. Normalized temporal growth rates obtained from spatial 
runs 2b (see Table 1) as a function of the number of grid points 
(in the longitudinal direction) used to describe one wavelength 
A x 2 of the linearly fastest growing mode. 



which increases exponentially in time during the linear devel- 
opment of the instability. 

For the spatial run, one cannot easily separate the spatial 
and time dependencies. For example, Fig. [7] shows an isocon- 
tour diagram of the mean transverse kinetic energy dU in the 
(?, x) plane. The range of the isocontour values is taken well 
above the amplitude of the initial perturbation and sufficiently 
below the amplitude of saturation of the instabilities, in or- 
der to address the linear regime. The ten different isocontour 
values that are plotted in the diagram are taken to be equally 
spaced between the minimum and the maximum value. Thus, 
the latter figure allows us to "catch" the growing instabilities 
as their mean transverse kinetic energy reach the level corre- 
sponding to the isocontour values. Each of the 1 1 bump-like 
structures shown in Fig. [7] indicates the presence of a grow- 
ing mode. Furthermore, the isocontours appear to become more 
and more clustered (in the (f, x) plane) as time goes on (recall 
that the isocontour values are linearly spaced). This is due to 
the exponential growth of the instabilities. In the same way, 
each "bump" in successive isocontours is shifted according to 
the phase velocity of the different modes. One can deduce the 
phase velocity of the different modes simply by measuring the 
slope of the characteristic curve joining the different isocontour 
curves. We estimated the phase velocity of the instabilities to 
be V ~ 0.65, which is close to the expected theoretical value 
of V p = Vq/2 (e.g. Drazin & Reid 1981). Moreover, a simple 
fitting of the transverse kinetic energy along the corresponding 
characteristic by an exponential function gives an estimate of 
the temporal growth rate. It is important to note that instabilities 
situated below the characteristic curve x = V p t (corresponding 
to a mode growing from x — at t — and shown as a dashed 
line in Fig. |7j are successively growing modes from the initial 
perturbation (situated at x = 0). On the other hand, the instabil- 
ities situated above this characteristic curve are the secondary 
instabilities which are driven by the primary modes of the in- 



stability train (as discussed in the previous section). As these 
secondary instabilities correspond to a non-linear driving phe- 
nomenon, we do not expect them to be correctly described by a 
linear theory. We deduce from Fig.|7]a temporal linear growth 
rate y w 0.125Vo/(2a) (see also Fig. 8). A very similar value 
can be deduced from our additional temporal run and also from 
Paper I. 

Furthermore, we checked the numerical convergence of our 
spatial results during the linear regime by repeating run 2 at 
different resolutions (see run 2b in Table 1). As we are only 
interested in the linear regime, we can safely use a smaller nu- 
merical domain to reduce the numerical cost. The measured 
temporal growth rates are plotted in Fig. [8] showing that us- 
ing 60 grid points per longitudinal wavelength is sufficient to 
achieve a good convergence. Moreover, note that taking only 
30 grid points is still acceptable, while allowing a considerable 
reduction of the numerical cost. Such "low resolution" run is 
used to investigate the long-term evolution (Sect. [33}. 

In conclusion, our study does not show any noticeable 
difference in the linear development of the KH instability in 
the spatial vs temporal approaches. We checked that this is 
also true during the early non-linear stages of our simulations, 
where a similar scenario with pairing/merging events between 
adjacent vortices leading subsequently to a large-scale disrup- 
tion occurs before x ~ 40/?/. However, this is no longer true for 
the further non-linear evolution of the system, as we will see in 
the following section. 

3.3. Stabilizing mechanism 
3.3.1. Description 

At later times in a temporal simulation, the jet configuration 
is able to reach a relaxed state which is a quasi-steady residual 
flow (Paper I). The final flow is thus drastically enlarged, invad- 
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20 40 60 80 



x 

Fig. 9. Axial velocity distributions from run 6 (same parameters as in run 2, but with a numerical domain twice as large) a long 
time after the large-scale disruption (f = 320). Dark regions correspond to low velocity values. A linear scale is used with values 
ranging from 0. to 1 .29. The unit of length is the jet diameter. 




Fig. 11. One dimensional cut at x — 3 of the axial magnetic 
field in Fig. [10] (plain line). The velocity profile is overplotted 
for comparison (dashed line). 



can see in the latter figure, slightly outside the jet's radius the 
magnetic field is amplified by a factor up to ~ 3.5 (i.e. ~ Ma/2) 
compared to its initial value. This is not true in the jet core, 
where the magnetic field has only slightly decreased. This new 
magnetic structure is observed to be maintained in time, and 
thus it changes the stability of the jet flow at later times. 

The question remains, what is the physical origin of this 
amplification? It has previously been shown that during the lin- 
ear and early non-linear evolution, the KH vortices (here of the 
instability train) are able to expel the magnetic field from the 
vortex center, stretching it and amplifying it around the vortex 
perimeter (e.g. Jones et al. 1997). A non-linear saturation then 
occurs when the magnetic field becomes locally dominant, i.e. 
when the field line tension is able to overcome the centrifugal 
force associated with the vortical motion. In our spatial simula- 
tions, we observe that the corresponding magnetic perturbation 
is able to travel backward up to the jet inlet, giving rise to the 
trail seen before. This is possible because the magnetic pertur- 
bation is traveling at the local Alfven speed in the jet's frame 
which has a local speed close to zero (see Fig. 1 1 at y « 1). 



ing the whole numerical domain, and most of the initial kinetic 
energy is transferred to the external medium. Consequently, the 
initial jet flow is not able to survive in temporal simulations. 

However, our spatial simulations indicate that a different 
scenario occurs at a later time. This is obvious in Fig. [9] (at 
t = 320), showing that the jet flow has been revived after 
the large-scale disruption. This figure was obtained from run 6 
when the dimension of the domain was twice as long compared 
to our standard run 2, in order to investigate the large-scale be- 
havior of the jet. In Sect. 13.11 we proposed the presence of a 
trail structure propagating in the upstream direction, as shown 
for example in the density map of Fig 3. 

This behavior correlates to a local strengthening of the 
magnetic field. This is clearly visible in Fig. [10] which shows 
a map of the magnitude of the magnetic field for run 2 at time 
t — 80 where the trail perturbation has just reached the jet inlet. 
In addition, Fig. QT] shows the corresponding transverse profile 
of the longitudinal magnetic field. The cut is made close to the 
inlet boundary (at x = 3), but a similar shape is also obtained 
along nearly the whole jet at this time (i.e. for x < 20). As one 



3.3.2. Influence of M A 

We investigated the influence of the choice of the initial mag- 
nitude of the magnetic field (measured using the Alfven Mach 
number Ma) on the final configuration. The results obtained for 
the final magnetic profiles are plotted in Fig.[T2]using the local 
Alfvenic Mach number defined as: 

M' A = V /V A (y) (5) 

where Vo is the velocity of the jet and VaOO = B x (y)/ ^p(y) is 
the local Alfvenic speed. This number is similar to the initial 
Ma, but uses the local value of the magnetic field intensity and 
plasma density. Figure [12] shows that a value close to 2 is lo- 
cally reached, irrespective of the initial Ma value. This means 
that the field is amplified up to a local maximum value corre- 
sponding to M l A = 2. Consequently, the amplification mecha- 
nism is rather robust as it works in a similar way for a rather 
large range of initial magnetic field amplitudes for transonic 
jets. 
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Fig. 12. One dimensional cut at x — 3 of the local Alfvenic 
Mach number M' A (see text) for different initial Alfvenic Mach 
numbers Ma = 5, 7, 10, 14 (from the bottom to the top) corre- 
sponding to the transonic runs 1, 2, 3 and 4 (see Tab. 1). 



3.3.3. Stability analysis of the final configuration 

A single shear flow layer embedded in a uniform magnetic field 
with an Alfven Mach number smaller than 2 is known to be 
fully linearly stable. However, our final 2D slab jet configu- 
ration is different (i.e. a double shear flow layer with a non- 
uniform magnetic field) and thus requires a stability analysis. 

We investigated the solutions of the dispersion relation de- 
rived by Hardee et al. for a magnetized slab jet (i.e. given by 
Eq. (1) in Hardee et al. 1992). An idealized equilibrium state 
is assumed, as the latter dispersion relation requires an equi- 
librium configuration with a vortex sheet (i.e. a discontinuity 
between the jet flow and the external medium). First, as an 
equilibrium state is not strictly reached at the end of our spatial 
simulation (see below), we calculated a neighboring equilib- 
rium having similar physical parameters. Second, we consider 
uniform magnetic fields inside the jet core (subscript "in") with 
A^A.in = 7 and in the external medium (subscript "ex") with 



Fig. 13. Temporal growth rate curves obtained by solving 
Hardee et al. (1992) dispersion relation with the vortex sheet 
approximation, for the initial equilibrium (upper curve) and for 
the final equilibrium (lower curve). The growth rates obtained 
using VAC code for 5 wavenumbers are also plotted (stars and 
diamonds are used for initial and final equilibria respectively). 



^A,ex = 2. The density ratio 77 = p m /p ex is taken to be unity 
and the sonic Mach numbers are consequently M s ; n = 1, and 
M s , ex = l.ll. 

The results are plotted in Fig. 13 for the temporal growth 
rate of the antisymmetric mode (which is the most unstable 
one). For comparison, the same stability curve for our initial 
uniformly magnetized configuration was also computed (in this 
case the parameters are 77 = pi n /p e x = 1, A^s.in = M sfiX = 1, and 
Ma.w = A^A.ex =7). One can see the enhanced stability effect 
for the final configuration (compared to the initial one) due to 
the enhanced magnetization of the external medium. However, 
a full stability is not achieved. As the latter stability study does 
not take into account the effect of the smooth transition veloc- 
ity profile at the jet interface (over a length scale a = 0.125), 
we additionally plotted the temporal growth rate of five wave- 
lengths, measured from additional temporal runs where a sin- 
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gle wavelength is driven unstable. The values obtained agree 
with the results deduced from the vortex sheet curves for wave- 
length A » a (or equivalent for small enough wavenumbers). 
However, small wavelengths (or equivalently ka < 0.7) become 
stable as expected from the effect of velocity gradients (Ferrari 
et al. 1982). The linear stability enhancement between the two 
equilibria is of the order of 50 percent, as measured when tak- 
ing the wavelength of the fastest linearly growing mode. 

Thus, according to our analysis, the linear stability of the 
final jet configuration is enhanced by a factor of approximately 
two due to the presence of a relatively strongly magnetized 
external medium (with a local Alfven Mach number of order 
2). However, this is only a partial stabilization, as unstable 
modes with non-negligible growth rates remain present. This 
is in agreement with the relatively short wavelengths seen at 
x m 20 in Fig. [9] Note that these short wavelength modes do 
not evolve into KH vortices, i.e. the non-linear evolution of the 
instability is different in the new configuration. 

As is the case for the residual large-scale oscillation of the 
whole jet with wavelength A ~ 30 (see Fig. 9 at x « 30), it is 
unlikely that this is due to the development of KH modes as 
our stability analysis predicts full stability for the correspond- 
ing wavenumber ka ~ 0.02. In fact, we have checked that, at 
the end of the simulation, an exact equilibrium is not yet fully 
reached and the system therefore oscillates when trying to relax 
towards it. 

3.4. Supersonic case 

We also performed a jet simulation having a supersonic flow 
with a Mach number M s = 3 (run 8 in Table 1) and show an 
axial velocity snapshot in Fig. [14] Two essential differences are 
evident when compared to the transonic regime. First, the jet 
perturbation is characterized by sinusoidal oscillations of the 
whole flow and not by vortex-like structures (observed in tran- 
sonic jets). This is due to the fact that the modes change charac- 
ter when Mf > 2 (in our case Mf = 2.75, Paper I). We deduce 
from the periodicity of the observed sinusoidal features that the 
wavelength of the fastest growing mode is A « 4 in our units 
(giving thus ka w 0.2), in agreement with previous linear sta- 
bility results (see lower curve in Fig. 2 of Paper I). In the latter 
stability study, the dominant mode is shown to be a KH surface 
mode. As concerns the growth rate, we deduce a value smaller 
than that obtained from the transonic run. This is not surpris- 
ing as the difference is mainly due to the fact that a perturbation 
convected by a supersonic flow is able to propagate further than 
by a transonic flow during the same given time. KH modes in a 
transonic flow (giving rise to two rows of vortices at the inter- 
face) translates into "surface" KH modes (in a supersonic flow) 
having longer wavelengths and lower spatial growth rates. In 
addition to the surface mode, internal modes become unstable 
in the supersonic regime (again when Mf > 2). These "body" 
modes are typical of two-shear and cylindrical layer configura- 
tions, as they become unstable by resonant reflection at the jet 
boundary (Ferrari 1998; Birkinshaw 1991). The latter modes 
are not expected to be dangerous for the jet survival. In our 
case, the dominant body modes are expected to have a wave- 



length A < 1 (in our units) and a rather low growth rate in 
comparison to the surface mode (see lower curve in Fig. 2 of 
Paper I), thus explaining why they are not visible in Fig. [14] 

The long-term evolution of our run shows that the previ- 
ously observed stabilization effect of the jet is now absent. This 
is not surprising, as the mechanism of magnetic field enhance- 
ment driven by the vortical deformation of the interface is not 
present for the supersonic flow. We also observe that the jet 
is disrupted in a similar way as obtained in temporal runs of 
Paper I, where shock dominated transients characterize the dis- 
ruption. Moreover, a propagation distance of order lOORj is ob- 
served in this supersonic case. 

4. Conclusion 

We can summarize our findings as follows. We have numeri- 
cally investigated the spatial development of the KH instability 
in jets embedded in uniform magnetic fields. We focus mainly 
on the transonic regime for weakly magnetized plasmas, for 
which the plasma-beta is much greater than unity. 

First, we confirm that a temporal approach provides a good 
approximation for modeling the initial stages of development 
of the KH modes. The latter assertion comes from a detailed 
comparison between our spatial runs and previous/additional 
temporal simulations. This is true for early times correspond- 
ing to the linear and early non-linear phases, ending up with a 
large-scale disruption of the flow. 

However, at later times the evolution differs drastically in 
our spatial vs temporal runs. Indeed, a temporal run is known to 
definitively terminate in a broadened and heated residual flow. 
As a consequence, a temporal transonic jet flow is not able to 
survive on an equivalent distance exceeding approximately 40 
times the jet radius. On the other hand, our transonic spatial run 
clearly shows that the initial jet flow is able to revive and main- 
tain its coherence at a later time after the large-scale disruption. 
A corresponding partial stabilization mechanism is identified. 
Indeed, the primary KH modes induce the formation of a sheath 
region of enhanced magnetic field situated around the jet core. 
This region is an envelope that primarily coincides with the 
region where KH vortices saturate, and progressively extends 
itself in the upstream direction towards the jet inlet. In this way, 
the stability of the whole flow structure appears to become sig- 
nificantly enhanced, as local Alfven Mach values between 2 
and 4 are reached in this region. This stabilization mechanism 
is not at work for supersonic flows (rolling-up vortices are ab- 
sent) which consequently do not survive KH instabilities. In the 
supersonic regime, the temporal and spatial approaches lead to 
the same conclusion: a rapid shock dominated disruption ter- 
minates the jet flow. 

Pioneering studies on spatial simulations of the develop- 
ment of KH instabilities in 2D unstable jets were performed 
by Hardee et al. (Hardee et al. 1992, 1995). However, these 
studies have focused on a different magnetic regime, as rather 
strongly magnetized configurations with a plasma-beta of order 
unity were considered. In a real jet (YSO or AGN for exam- 
ple), the plasma parameters are not known even if one expects 
an average plasma-beta of order unity. For example, it cannot 
be excluded that a high value of the plasma-beta of order 100 is 
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Fig. 14. Axial velocity distribution in the supersonic case (run 8 : M s = 3, Ma = 7) at t = 300. The unit of length is the jet 
diameter. 



reached in a central region of the jet core (see Fig. 6 in Lery et 
al. 2000). Thus, before reaching the asymptotic region where 
the flow is supersonic, the plasma of a real jet is accelerated 
in an intermediate transonic region. The partial stabilization 
mechanism presented in this study could thus help to explain 
the jet survival in the transonic regime. This is possible if the 
transonic phase lasts long enough for the stabilization effect to 
operate. In the opposite case, the flow is anyway stable as the 
KH modes do not have enough time to grow. For example in the 
accretion/ejection model of Casse & Keppens (2002, 2004), ve- 
locity curves show that the jet is accelerated into the supersonic 
regime on a very short-length scale, of the order 50 times the 
inner radius of the accretion disc. As the jet in their simulation 
has a radius of approximatively 20 inner radius, it shows that 
according to this model the transonic phase has a very short 
duration. Unfortunately, our stabilization mechanism does not 
work in the supersonic regime, for which the problem of the 
jet survival remains unsolved in spite of many recent attempts 
(Paper I; Hardee & Rosen 2002; review by Hardee 2004). 

Note also that in the present study, we use idealized equi- 
librium profiles for the jet configuration to facilitate the inter- 
pretation of the results. It would be interesting to test the ro- 
bustness of the stabilizing mechanism in the presence of more 
realistic profiles. For example, the velocity profile expected in 
astrophysical jets is probably much more complex, consisting 
of different components (e.g. Frank et al. 2000). This latter in- 
teresting question is left to a future investigation. 

Finally, it will be of interest to study whether the results ob- 
tained in this study apply also in the 3D case. We do not expect 
strong differences, mainly because there is a certain correspon- 
dence in the mode nomenclature between 2D and 3D config- 
urations. The only known difference is the enhanced turbulent 
aspect in 3D due to the presence of typical 3D modes and non- 
linear couplings (Paper I). 
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